Stable and unstable vector dark solitons of coupled nonlinear Schrodinger equations. 
Application to two-component Bose-Einstein condensates 
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Dynamics of vector dark solitons in two-component Bose-Einstein condensates is studied within 
the framework of the coupled one-dimensional nonlinear Schrodinger (NLS) equations. We consider 
the small amplitude limit in which the coupled NLS equations are reduced to the coupled Korteweg- 
de Vries (KdV) equations. For a specific choice of the parameters the obtained coupled KdV 
equations are exactly integrable. We find that there exist two branches of (slow and fast) dark 
solitons corresponding to the two branches of the sound waves. Slow solitons, corresponding to the 
lower branch of the acoustic wave appear to be unstable and transform during the evolution into the 
stable fast solitons (corresponding to the upper branch of the dispersion law). Vector dark solitons 
of arbitrary depths are studied numerically. It is shown that effectively different parabolic traps, 
to which the two components are subjected, cause instability of the solitons leading to splitting of 
their components and subsequent decay. Simple phenomenological theory, describing oscillations of 
vector dark solitons in a magnetic trap is proposed. 
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I. INTRODUCTION 



It is well established that in a one-component Bose- 
Einstein condensate (BEC) with a positive scattering 
length, which has cigar-shaped geometry, one can gen- 
erate dark solitons Experimental generation of two- 
component BEC's of different hyperfme states of rubid- 
ium atoms in a magnetic trap Q and of sodium atoms in 
an optical trap stimulated theoretical studies devoted 
to the meanfield dynamics of multicomponent conden- 
sates. As in the one-component case special attention 
was devoted to existence of solitary waves in such sys- 
tems. When a condensate is cigar-shaped and has rela- 
tively low density, i.e. when the healing lengths of the 
components are much larger than the transverse dimen- 
sion of the condensate and much less than its longitudinal 
dimension, the transverse atomic distribution is well ap- 
proximated by the Gaussian ground state and the system 
of coupled Gross-Pitaevskii (GP) equations, describing 
the mixture (see e.g. 0), can be reduced to the coupled 
one-dimensional (ID) nonlinear Schrodinger (NLS) equa- 
tions (see e.g. [BJ). The respective models were a subject 
of recent theoretical studies. In particular, coupled large- 
amplitude dark-bright solitons have been reported in km ; 
bound dark solitons have been numerically studied in pj , 
where it has been found that creation of slowly moving 
objects is possible; a diversity of other bound states has 
been generated numerically in 

The present paper aims further analytical and numeri- 



cal study of dark solitons in two-component BEC's. Main 
distinctions of the situation considered here compared to 
the previous research are as follows, (i) We consider vec- 
tor dark solitons, i.e. states where two components move 
with equal or approximately equal velocities, (ii) We do 
not impose the condition of equality of the nonlinear coef- 
ficients, as a necessary condition, allowing one to reduce 
the problem to the exactly integrable one - to the so- 
called Manakov problem, for which vector dark solitons 
are known [jj. (iii) In the small amplitude limit we pro- 
vide analytical description of the phenomenon reducing 
a system of coupled NLS equations to a system of cou- 
pled Korteweg-de Vries (KdV) equations, what allows us 
to predict existence of two types of vector dark solitons, 
moving with different velocities, (iv) Finally, we study in 
detail the effect of the magnetic trap on the dark soliton 
dynamics. We show that, due to difference of its effect 
on different component, a magnetic trap leads to splitting 
of the components and subsequent destruction of vector 
dark solitons. 



II. STATEMENT OF THE PROBLEM AND 
PHYSICAL PARAMETERS 

Evolution of a two-component BEC composed of dif- 
ferent hyperfme states is described by the coupled GP 
equations (j — 1,2) Q 



dt 



2m rn ^-^ 
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where Vj(r) — —uj j (A x + rj_), V&j are the macroscopic 
2 

wave functions of the states, scattering lengths of 



2 



the respective interactions - it will be assumed that they 
are positive, luj are transverse linear oscillator frequen- 
cies of the components, and A is the aspect ratio of the 
condensate. Respectively, Afj — J \ 2 dr is the number 
of atoms of j-th component and JV = N\ +N2 is the total 
number of atoms. 

In the case of an elongated trap, when A is small 
enough, and when densities of the both components are 
low enough, one can employ the multiple-scale expansion 
method in order to reduce the original 3D system Q to 
the homogeneous coupled ID GP equations (for details 
of derivation see e.g. y|) 



X*2 



Xi|$i| 2 $i +x|$ 2 | 2 *i, 

X|$l| 2 $2 +X2|* 2 | 2 *2- 



(2) 



III. SMALL AMPLITUDE DARK SOLITONS 
A. Rescaled system of equations 

An initial value problem for system does not al- 
low solution in a general case, except in the special limit, 
which is known as the Manakov system and which is dis- 
cussed below [see (1121 the respective discussion]. Some 
information about possible solutions, is however avail- 
able in the small amplitude limit, where the coupled NLS 
equations are reduced to the coupled KdV equations. 

In this approximation a dark soliton evolves against a 
background what means that (J2J is considered subject to 
the boundary conditions 



Here 



lim \<f 3 \ 2 = P 2 



(5) 



Xi 
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X2 = r- 
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$j, T and X are the dimensionless wave function en- 



velops, slow time and slow coordinate, respectively, u> 

9ij = 47r "^ a ' 3 , and aj = y 
eter of the problem is defined as 



h 



where p 2 are properly normalized densities of the com- 
ponents. Taking this into account as well as a large num- 
ber of free parameters of the problem, is convenient to 
scale out the boundary conditions for the sake of per- 
forming the small-amplitude reduction mentioned above. 
To this end we introduce the total dimensionless den- 
sity p 2 = p\ + P2, and rescale the variables as follows: 
tpj = t = XP 2 T, and x = ^/xpX. Then Eq. © is 



Pi 



The small param- rewritten in the form 



S = 



jTrMaiiA 1 / 2 
ai 



(4) 



Besides asymmetry of the trap and weakness of the 
two-body interactions, expressed by the smallncss of A 
and 5, respectively, our model assumes equality of the 
aspect ratios of the components (A does not depend on 
j). Meantime we emphasize that the linear oscillator fre- 
quencies u>j may be very different (say, in the experiment 
reported in Ref. ^(| , the relation between the frequencies 
was lo — i/2). As a result, even for initially equal scat- 
tering lengths, the effective nonlinearities \ji defined in 
Eq. JjU, become different because of different transverse 
distributions of the components. 

To estimate a typical value of the parameter 8 we con- 
sider a binary condensate of two hyperfine states of ru- 
bidium atoms in a trap with the mean value of the trans- 
verse oscillator frequency 2ir x 200 Hz and the aspect ratio 
A = 10~ 4 (this corresponds to the 1 pm and 100 pm of 
the transverse and longitudinal linear oscillator lengths) . 
Taking an w 1 nm (here we take into account that the 
scattering length can be varied by using Feshbach res- 
onance) and assuming the mean atomic density to be 
n ps 10 12 cm~ 3 we obtain 5 » 0.05. We also point out 
here that the respective healing length £ (for numerical 
estimates we assume that the healing lengths of the both 
components are approximately equal) is estimated to be 
of order of 4 pm. 



idtipi = -d 2 V>i + (t/i|V'i| 2 + 



cos 2 a\ip 2 \ 2 )ipi 



-^2 + (sin^a|^ 1 | 2 + C/ 2 |V'2| 2 )V'2 



(6) 



where Uj — 



cos a 



= sin a = and thus the 



'3 ~ xp' 2 

parameter a determines the relation between the unper- 
turbed densities of the components: tana = The 
boundary conditions now acquire the desirable form 



lim 



1. 



(7) 



B. Sound propagation 

Small amplitude dark solitons, which from the phys- 
ical point of view represent packets of acoustic waves, 
for which weak nonlinearity and weak dispersion are bal- 
anced, propagate against a background with a speed close 
to the group velocity of the sound (see e.g. as 
well as consideration below). The background in our case 
is computed from © and J7J to be ipj = exp(—i£jt) 
where 



£ 1 = U\ + cos 2 a, £2 = U2 



sin 2 a . 



(8) 



Designating frequency and wave vector of a sound wave 
as n and K, respectively, i.e. considering a solu- 



tion of © in a form tpj = 1 
Cj exp [—i(Kx — fli)], where bj 
stants, I61 



-I- bjesxp[i(Kx — Clt)] + 
and Cj are small con- 
<C 1, one finds the two branches (upper 
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with sign "+" and lower with sign "— ") of the spectrum 
of the acoustic waves 



0± = K\Jk 2 + U 1 + U 2 ± y / (L/ 2 -C/ 1 ) 2 +sin 2 (2Q). (9) 

Since dark solitons will be constructed against a static 
background we are interested in the limit of long wave- 
lengths, where the group velocities are given by: 

,. dn± 



= \ju 2 + U 1 ±^{U 2 - t/i) 2 + sin 2 (2a). (10) 

It follows from @, that the lower branch is stable subject 
to the condition 

A = UiU 2 - sin 2 a cos 2 a > (11) 

and thus, the consideration below will be restricted only 
to this case. We notice that this constrain corresponds 
to the condition of the thermodynamic stability of the 
condensate with the only difference, that it is written for 
the effective nonlinearities (rather than in terms of the 
coefficients gy, see e.g. Q). 

An interesting feature, relevant to the next considera- 
tion, is that the group velocity of the long-length excita- 
tions of the lowest branch, i.e. i>_, becomes zero at 

sin 2 (2a) = AUiU 2 . (12) 

Then the matrix of effective nonlinearities A is degener- 
ated. This situation corresponds to the Manakov system, 
which is an integrable limit of the coupled NLS equations. 



could provide the inequality 8 <C e 5 <C 1, what is not fea- 
sible for real experimental situation, where typical values 
of S are of order of 0.1 -j- 0.01, as it has been mentioned 
above]. Thus the present section although being inter- 
ested from the viewpoint of the small-amplitude limit 
of the coupled NLS equations, cannot be directly inter- 
preted as a theory of small amplitude dark BEC solitons 
(the issue which remains to be an open problem). The 
significance of the results obtained below for the mean- 
field theory of the BEC is in the indication on new types 
of the solutions which cannot be directly obtained from 
the system © (or ©). 

There exists one more limitation for the practical use 
of the small amplitude limit. It is related to the fact that 
small amplitude solitons are wide and their width may 
be comparable with the longitudinal extension of the con- 
densate. Indeed, the characteristic size of a (non-small- 
amplitude) dark soliton is of order of a healing length £j 
(as above healing lengths of the both components £1,2 are 
considered of the same order). This means that a small 
amplitude soliton has a width of order of £j/e ^> £j. 
The small amplitude expansion fails at the boundaries 
of the atomic cloud, as the density approaches zero in 
those domains. Thus, for the validity of the theory one 
must require the longitudinal size of the condensate (i.e. 
aj/^f\) to be much bigger than the width of the soli- 
ton aj/VA ^> £,j/e or in other words one must impose 
a condition e 3> V\£,j/aj. The obtained constrain is , 
however, not as strong as the previous one. In particu- 
lar, the estimates provided the end of the Sec. II give 
now e > 0.04. 



D. Coupled KdV equations 



C. A comment on the small parameter 

Let us now turn to the analysis of small amplitude dark 
solitons. To this end we employ the idea due to Ref. ^l| 
about possibility of the multiple-scale reductions between 
NLS and KdV equation, noticing that the KdV equation 
can be obtained as a small amplitude limit of the NLS like 
equation also in nonintegrable limit and with arbitrary 
(intensity dependent) nonlinearity |T^. 

Before going into details of calculations we make the 
following observation. Derivation of the dynamical equa- 
tion for small amplitude waves is based on introduction 
of the small parameter of the problem, which we will be 
designated as e [see l|Ti]l. (fTIfl) below]. Then the resulting 
evolution equation (see (|18(l below) appears in the order 
e 5 . If we now recall that the ID reduction of the coupled 
GP equations is mathematically justified by smallness 
of 8, introduced in (0J, and that the respective ID NLS 
equations (0) appear in the S 3 order, while terms of order 
of 8 4 are neglected, we conclude that strictly speaking the 
resulting KdV limit is not applicable for the description 
of low-dimensional BEC's [it would be valid only if one 



We look for a solution of system © in a form of a 
small amplitude excitation of the background exp(— i£jt), 
which moves with a velocity close to one of the speeds 
v± given by (|1C)|> . The respective analytical ansatz reads 

tpj = Qj(C,T)exp(-i£jt + i(pj(C,r)) (13) 

where the amplitude, Qj(£,r), and the phase, ipj(Ct T )i 
are represented in forms of the expansions 

Qj((,T) = l + £ 2 9 ,(C,T)+eV(C,T) + ... , (14) 
<pj (C, t) = (C, r) + e 3 ^-! (£ r) + • • ■ (15) 

and the new slow variables are given by 

C = s(x — vt), 



([/i + [/ 2 )[2(£/i + [/ 2 )v 2 -4A] ' 

Hereafter v is either v+ or u_, depending on the branch 
under consideration. One verifies, that subject to con- 
dition 1|11[1 the denominator in the definition of the slow 
time r, which is introduced for the sake of convenience, 
is always positive. 
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Our aim now is to derive evolution equations for 
t) and <pj{C, t), which will describe evolution of the 
small amplitude dark solitons. Respectively, we impose 
the boundary conditions 



lim q j ((,T) = 0, lim <^-(C, r) = (^±, (16) 



with cf>j± being constants, and consider equations of dif- 
ferent orders of e. 

While equations of the zero and the first orders are 
satisfied identically, it follows from the equations of the 
2-nd and 3-rd orders with respect to e (see Appendix A 
for the details), that the amplitudes and the phases of 
the excitations of the two components are linked by the 
relations (j = 1, 2) 



q 3 = ~d£(l>j. 



(17) 



The obtained formula reveals essential difference between 
integrable and nonintegrable versions of the coupled NLS 
equations. In the former case a zero value of the group 
velocity of the lowest branch does not allow existence 
of small amplitude solitary pulses. In other words in 
the integrable case there exists only one branch of dark 
solitons. 

The condition of compatibility of the equations of 4-th 
and 5-th orders in the small parameter e results in the 
coupled KdV equations (see Appendix A): 



d T q 3 

where j, k, I = 1, 2 



daflkqi 



d 3 c (3^q k = 



(18) 



7 X U = Qv 2 Uf + 8v 2 U 1 U 2 + 4w 2 A - 12Z7 X A , 

7 f = 2 cos 2 a [v 2 {2Ui + U 2 ) - 4A] , 

7i 2 = if = cos 2 a [v 2 {Ux + 2U 2 - sin 2 a) - 2 A] , 

7 2 n = 2 sin 2 a [v 2 ^ + 2U 2 ) - 4A] , 

7 f = 6v 2 Ui + 8v 2 LW 2 + 4w 2 A - 12J7 2 A , 

if = if = sin2 ot [v 2 (2Ui + U 2 - cos 2 a) - 2 A] 



are the effective nonlinearities, 



Pl=vA-— {U 1 +2U 2 ), (5\ 



3 



cos a . 



Pi =vA-- (2U, + U 2 ), (3\ = - sin" a 

arc the effective dispersions, and the Einstein summation 
rule over repeated indexes is used. 

Let us return to discussion of the integrable limit (|12|l , 
where according to the discussion following (|17|) . one 
must take v = i> + = ^/2(Ui + U 2 ) (i.e. the upper branch 
only) and A = 0. In this case, the coefficients of the 



coupled KdV equations become 7? = {U\ + U^-fj* and 



Pi = ^{Ui + U 2 )Pl where 



(3l = -2(£7i + 2£7 2 ), 0\ = 1 + ^1 - AUxU 2 , 
/3 2 = -2(2C/ 1 + C/ 2 ), p\ = 1 - y/1 - 4U X U 2 , 

7l 11 = i(2/3 2 -/3 1 1 )(5/3 1 1 + 2/3 2 2 ), 
7i 12 =7i 21 --k 2 (/3i 1 +/3 2 1 ), 



if 



PlPf 7" = , 



1 5i 



7 2 " = 7 21 



1 - 



7 2 2 = -(2/3 1 1 -/3 2 )(5/3 2 + 2/3 1 1 ) 

(notice that now U\U 2 < 1/4). The respective system of 
the coupled KdV equations is integrable. 



E. Small amplitude dark solutions 

A specific particular solution of Eq. I|18|) can be 
searched in the form 

V 



qi = q-2 



cosh 



f«(C-«^) 



(19) 



where w, and r\ are real constants to be determined and n 
is a constant parameterizing the problem. Substitution of 
ansatz <|19fl in Eq. I|18|) gives the condition of the equality 
of the chemical potentials 



£1 =£2 



(20) 



Next, the parameter 77, characterizing the width of the 
soliton, is computed to be 



2k 2 



V = 



1 + Ui + U 2 



(21) 



The parameter w, which characterizes the soliton velocity 
in the frame moving with the speed of the sound, is given 
by 

w = w + = -4([/ 1 + C/ 2 )(1 + [/i + [/ 2 )k 2 (22) 
for the upper branch and 

w — w- = -4(f/i + U 2 



1) 2 K 2 



(23) 



for the lower branch of the spectrum. Both w± are neg- 
ative, what means that the solitons move with veloci- 
ties smaller than the sound velocities of the respective 
branches. 

As it follows directly from ansatz (|19fl . solutions de- 
scribing different branches correspond to equal distribu- 
tions of atomic densities. Meantime they correspond to 
different phase differences 

A<pj = <t> H - (24) 

at the infinities [see Eqs. i|17|) an d <|lt>[l ]: the lower- 
branch soliton "separates" domains with larger difference 
of chemical potentials. 
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F. Numerical results 

As it has been mentioned above, although solutions of 
the coupled KdV equations represent a good approxima- 
tion for the exact solutions of the coupled NLS equation 
in the small amplitude limit, strictly speaking they can- 
not be considered as satisfactory, when applied to the 
dynamics of a two-component BEC in an elongated trap. 
Meantime, the obtained vector soliton (|19fl can be em- 
ployed as an initial condition for numerical generation 
of dark vector solitons of the coupled NLS equations 
having small but finite amplitudes. Such numerical study 
is performed in the present section. 

Taking into account the existence of two kinds of small 
amplitude solitons we address the question about persis- 
tence of the excitations at finite, but still small, ampli- 
tudes, as well as their stability withing the framework of 
system (jSJ. 

An exact vector dark soliton corresponding to the up- 
per branch of the linear spectrum reads 



where 



i v o + y 1 — t arm [ K ( x ~ Vot)] 



,-iSjt 



4k 2 



"0 



Vo 



(25) 



(26) 



and relation l|20[). which also can be rewritten as 
2 cos 2 a = 1 — Ui + U2, is taken into account. By ex- 
panding H25fl in the Taylor series in terms of the small 
parameter n one verifies that the leading orders trans- 
form into the upper-branch dark soliton described by the 
formulas (T3J, iJTSJ, (jHJ, (|2TJ>, and J22J), where the for- 
mal small parameter e is substituted by one and k is 
interpreted as the small parameter of the problem. 

In Fig^a) the trajectories of the centers of vector dark 
solitons are shown for three different initial conditions: 
the exact dark soliton (|25l) (line 1), the approximate dis- 
tribution l|19fl corresponding to the upper (line 2), and 
lower (line 3) branches. The centers of the solutions are 
defined as coordinates of the absolute minima of \ipj\ 2 '- 

they are designated as x^ in and x^in f° r the exac t and 
two approximate solutions, respectively. 

Exact solution H25fl appears to be stable, and its ap- 
proximate counterpart given by l|19f) undergoes small (in- 
visible on the scale of Fig^l deformation. In the numer- 
ical simulation with the "approximate" initial condition 
corresponding to l|19|) the effective small parameter n is 
0.1. This gives the amplitude difference with the exact 
solution to be of order of k 4 = 10~ 4 what explains differ- 
ence between the lines 1 and 2 in Fig^b). The solution 
corresponding to the lower branch of the spectrum starts 
to move with the velocity w_, but during initial interval 
of time changes significantly: it decays into two local- 
ized pulses moving in opposite directions, as it is shown 
in Fig[3 The forward moving part represents the "un- 
stable" dark soliton. Its amplitude continues to change 



|v|>fc,,)l 




FIG. 1: (a) Trajectories of dark solitons. Lines 1, 2 and 3 



correspond to dynamics of x^ in , x^ n and x y m - n computed 
numerically. The dotted line corresponds to the theoretical 
prediction of dynamics of the x ~^ n . Theoretically predicted 
values of x^ in and x^In practically coincide, what makes 
them indistinguishable in the figure. Parameters are Ui = 2, 
U2 = 1.5, a — 7r/3, k = 0.1, e = 1. (b) Dynamics of the 
minima of densities of the dark vector solitons: the grows of 
the minimum shown by the line 3 corresponds to shallowing 
the soliton solution. 



,(+) 



„(-) 



during the evolution (see Fignjb)) and its velocity ap- 
proaches the group velocity of the upper (stable) branch 
(see FiglUa)). 

The observed instability of the lower branch dark soli- 
ton is especially interesting in view of the recent results 
reported in [13J, where it was found that in the 3D cou- 
pled GP equations, describing spinor condensates, there 
exist subsonic (i.e. having velocities below the sound 
speed of the lowest branch) solitary wave complexes, 
most of which are suggested to be stable. In our case 
a subsonic vector dark soliton is unstable, and the sta- 
ble one has the velocity between the two sound speeds: 
V- < Vo < v + . 

In order to understand this phenomenon, let us con- 
sider the energy of the system 



E — Ei + E2 + Ei n t 
which we write down in dimcnsionlcss variables 



Ej = J dx 



^(l^l 2 -!) 2 



(27) 



(28) 



is the energy of the excitation of j-th component and 

£U = i / dx(|^| 2 -l)(|^ 2 | 2 -l) (29) 

is the energy of interaction of the components. Next, we 
recall that for a given phase differences between the in- 
finities, Aipj [see definition Q24p]. there exist two types 
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FIG. 2: Profile of the unstable lower-branch soliton at t = 
(thin line) and at t = 10 (thick line). Parameters are the 
same as in Fig0 



of dark solitons, corresponding to two branches of the 
group velocities. Let us assume now that the respective 
lower and upper branch solitons are characterized by pa- 
rameters k_ and k+, respectively. Let also initially the 
low branch soliton is excited. Since it does not represent 
an exact (but only an approximate) solution it starts to 
deform with time. Such a deformation has a constrain: 
the phase differences Aipj are preserved. In the case of 
a small amplitude soliton, given by {15} . one computes 
from Eq. JT7J and definition lf24"|) 



2 3/2 



1 + Ui + U 2 



3/2 
K+V + = 



2 3/2 



1 + Ux + U 2 



k-v 3 / 2 . (30) 



On the other hand, one can compute the energy in 
terms of the parameters of the solution by substitution 
of (Uni) and |2T} into expressions (|T7|l - (t2T?|) : 

_ 164 [4 + 2(1 + Ui + u 2 )] 



From (|30|) it follows that k + v 
V— < v+ one finds that n + < k 



3/2 



3/2 c . 

. Since 



and hence E+ < E-. 
Thus what we observed in the numerical simulations is 
the transformation of a high energy vector soliton in a low 
energy vector soliton, accompanied by quasilinear modes. 



IV. EFFECT OF THE PARABOLIC TRAP 

Let us now turn to the situation when the trap is not 
long enough in the axial direction, such that reflections 
of a soliton from potential walls can happen. In this case 
the trap potential must be included explicitly in the equa- 
tions (instead of including it into the normalized ground 
states, see e.g. 14]). This leads to coupled ID GP equa- 
tions © as follows 



idtipi = -flgVi + { V W + Ui\^i\ 2 + cos 2 a|V>2| 2 )V>i 
id t f 2 = -d^ 2 + (yW + sin 2 a|Vi| 2 + U 2 \^ 2 \ 2 )^ 2 



Here v\ — X/xP 2 an d v 2 = wv\ are effective strengths 
of parabolic traps. Without restriction of generality in 
what follows we consider the situation where v± <v 2 . 

In the case of equal effective frequencies, v\ — v 2l and 
subject to the specific initial condition (-01 =const- - 02) J 
effectively reducing the system to a single NLS equa- 
tion, a vector dark soliton oscillates with the frequency 
V2 times smaller than the frequency of the parabolic 
trap 0, 0| . If however strengths of the parabolic po- 
tentials of the two components are different, behavior of 
the soliton changes dramatically. There exist two com- 
peting factors, which determine the dynamics. On the 
one hand, each component, affected by its own parabolic 
potential, "attempts" to oscillate with its own frequency. 
On the other hand, attractive interaction between the 
components, described by functional H29J) . forces them to 
oscillate with the same frequency. 

The dynamics emerging from the competition of the 
factors mentioned above is shown in Fig|3J In the case of 
sufficiently wide potentials and small difference between 
their strengths (FigJSfa)) the first and second compo- 
nents (thin and thick lines correspondingly) initially os- 
cillate with approximately equal frequencies. 
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FIG. 3: Trajectories of the components of the vector soliton, 
x m in,j (j = 1,2), corresponding to the first (thin line) and 
to the second (thick line) components for parameters ft, v\ 
and vi being respectively: (a) 1.05, 0.15, and 0.2; (b) 1.05, 
0.1, and 0.2; and (c) 1., 0.1, and 0.2; (d) 1.05, 0.2, and 0.4. 
Dotted lines show trajectories of the components in the re- 
spective traps when interaction between the components is 
absent. The other parameters are U\ — 2, U2 = 1.5, and 
a — 7r/3. 



The dynamics observed in the case when initial soliton 
velocities are fixed and difference between strengths of 
the traps is increased is shown in FiglSfb). One observes 
more significant separation of the components. The dy- 
namics, however still resembles periodic motion. After 
subsequent increase of the initial kinetic energy of the 
soliton, by means of increasing its initial velocity, more 
visible splitting is observed, as it is shown in FigUJc). Af- 
ter the first period the soliton splits and each component 
starts to oscillate with its own frequency. 

By increasing the frequencies of the magnetic trap 
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but keeping constant the ratio between them (this cor- 
responds to passage from Fig^b) to FigEfd)) one ob- 
serves fast splitting the components whose trajectories 
show rather independent behavior. The second compo- 
nent, which exists in effectively more narrow trap, does 
not display periodic motion. 

To understand qualitatively the described behavior, we 
introduce coordinates of the centers of mass of the com- 
ponents (j = 1, 2) 



explicitly. It appears to be a function of X- only (i.e. 
independent on X + ). 

Let us now take into account that \X-\ <c + | [see 
Figs. |3| and 01 panels (a) and (b)]. Then, in the leading 
order the terms with X_ can be neglected in (|34[1 result- 
ing in a simple equation for the center of mass of the 
condensate: 



X 4 



nix. 



0. 



(36) 



x\ipj(x)\ 2 dx, Nj 



\tpj{x)\ 2 dx. (33) 



We also define the coordinate of the center of mass of 
the whole condensate X + — (NiXi + N 2 X 2 )/N, where 
N = N% + N 2 , and the distance between the two centers 
of masses: X_ = X% — X\. We emphasize that, strictly 
speaking, X\ 2 do not describe trajectories of the dark 
solitons (see the discussion in 0). One however could 
expect, that when the vector soliton splitting is small 
enough (i.e. when X- <C 1) the relations among the 
frequencies of the two-component problem are approxi- 
mately the same as in the one-component case. That is 
why we concentrate on dynamics of X±. 

Differentiating Xj with respect to time and using © 
we compute 



+ ^ {N lV l + N 2 v 2 2 ) X+ 



4 
N 



N 2 \ 2 A7iA7 2 2 



N 



— (sin 2 a — cos 2 a) 
N v ; 



,d\ip 2 



X. 



-dx, 



X. 



N2\ 2 ,N 23 
N J 2 N 2 



dx 

X- -v\)X 



(34) 



= 2 



sin 2 a cos 2 a 



Nn 



dx 



-dx . 



(35) 



Here 



tf + = -(N lV 2 + N 2V 2 ) 



(37) 



In other words Q s = Q + /\/2 is expected to be the main 
frequency of the oscillations of the vector soliton, where 
we introduce the factor y/2 conjecturing that the relation 
between the frequencies of the condensate and the soli- 
ton is the same as in the one-component case. To com- 
pare this estimate with the direct numerical simulations 
shown in FigOJa), we use the numerical values N\ 12 
and N 2 w 7.1 and obtain from (|37|l : the frequency 
ft s m 0.24 and the period of oscillations T s « 26.1. The 
numerical value of the period subtracted from Fig|3fa) 
is T n um,s ~ 26.3 what is in a remarkable (taking into 
account the character of the approximations) agreement 
with the theoretical prediction T s . 

In order to analyse the characteristics of the splitting 
of the soliton, we hold the above assumptions and rewrite 
Eq. I)35|l in the form 



X_ + Q 2 X_ 



dUeff{X.) 

dx_ 



4 {v\ - v 2 2 ) X + (38) 



where the effective potential is given by 



. 2 

sm a 



N 2 
sinh X- 



X_ coshX^ 



smb. 6 X. 



(39) 



Hereafter overdots stand for derivatives with respect to 
time. 

While equations 134|) and l|35l) are exact, they are not 
closed, and in order to make use of them we have to make 
some approximation. Below, in Figs. Efa),(b) we show 
that, if the difference v 2 — V\ is small enough, the shapes 
of the components are preserved for relatively long tem- 
poral intervals, even when the splitting is not negligible. 
Basing on this observation we assume that the compo- 
nents of the vector soliton preserve their shapes and only 
change their velocities. This allows to us describe each 
component in the vicinity of the soliton by formula 125|) 
where V$t is substituted by Xj(t). We also restrict the 
consideration (this time for the sake of simplicity only) 
to zero initial velocities: vq = 0. Then the integral in 
the right hand sides of l|34|l and l|35|l can be computed 



U e ff(X_) describes additional confinement of the rela- 
tive motion of the components due to the attractive in- 
teraction between the components, which is imposed si- 
multaneously with the parabolic trap characterized by 
the frequency 



n 2 _ = 4 



1 - 



N2 

N 



N- 



2 ,.2 



N 



(40) 



Eq. (|38J) is nothing but a nonlinear oscillator in a trap 
made up of the parabolic potential fi 2 _X 2 /2 and of the 
nonlinear potential U e tf(X—), driven by the periodic 
force, 4 (y\ — i/f) X + (t), originated by the oscillation of 
the condensate as a whole. 

To simplify the next consideration we take into account 
that X- is relatively small [see FigQJa)] and substitute 
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dU e ff(X_)/dX_ by the first term of its Taylor expan- 
sion. This results in a modified linear frequency, which 
is determined as 



j Q 2 _ 



^£/(0) 



(41) 



Using the data of FigEJa) we compute fi_ w 0.59, 
what corresponds to the period of modulation T_ = 
10.65, while the numerical value subtracted from the 
figure is T„„ m ._ 11.8. Thus we again observe good 
agreement between the simple theoretical estimates and 
the numerical results, what shows that the simple model 
given by Eqs. and provides adequate descrip- 
tion of the vector soliton dynamics in a parabolic trap 
whenever the splitting between the components is small. 



f t 


A (a) 
i i- 


1,1, 






FIG. 4: Dynamics of x m i n ,i — x m in,2 corresponding to FigEl 
The panels (a)-(d) correspond to the respective panels in 
FigH 

Interaction between the components of the vector soli- 
ton and soliton interaction with the confining potential in 
a general case lead to deformations of the profiles of the 
soliton components, what is shown in Fig|3]for the times 
tf corresponding to the final time of the evolution de- 
picted in Fig|3| Already after few oscillations the vector 
soliton decays, except in the case when both components 
oscillate in effective traps with equal strengths (i.e. when 
the vector soliton behaves like a single-component dark 
soliton 



V. CONCLUSION 

To conclude, we have investigated dynamics of vector 
dark solitons governed by one-dimensional coupled non- 
linear Schrodinger equations. In the homogeneous case 
and in the small amplitude limit, when the vector soliton 
propagates with the velocity close to the speeds of the 
sound, a stable vector soliton has a velocity close to the 
higher velocity of the sound and exceeding speed of the 
slow sound. The respective subsonic (i.e. having a veloc- 
ity lower than the speed of the slow sound) dark vector 
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FIG. 5: Density profiles |V>j| 2 of the first (solid lines) and 
second (dashed lines) components at time tf corresponding 
to FiglHl In (a)-(c) tf = 40 and in (d) iy=20. 



soliton is unstable. The both cases are described by the 
coupled Korteweg-de Vries equations. When the group 
velocity of the lower branch of the sound dispersion re- 
lation becomes zero, what corresponds to the integrable 
Manakov system, the coupled Korteweg-de Vries equa- 
tions obtained in the small amplitude limit are integrable 
(to the best of authors knowledge, such system has not 
been reported in the literature, so far). 

Including a parabolic trap into consideration changes 
behavior of vector solitons dramatically, leading in a gen- 
eral case to their decay, what is explained by different 
eigenfrequencies of the two components. If meantime 
the effective traps for the both components have close 
frequencies, during initial times the dynamics of the vec- 
tor soliton can be qualitatively (and also quantitatively, 
with rather good accuracy) described by the oscillatory 
motion of the soliton center with the frequency given by 
l|37p. Relative dynamics of the components, when split- 
ting is small enough, can be interpreted as an oscilla- 
tor driven by a periodic force. Large difference between 
strengths of parabolic traps or large initial soliton ve- 
locities cause instability of dark vector solitons leading 
to their splitting and subsequent decay, thus preventing 
possibility of the long-time dynamics. 
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APPENDIX A: DERIVATION OF THE COUPLED 
KDV EQUATIONS 

Substituting (JT^J, iJTSJl and (HJ in ©, and gath- 

ering terms of the same order of the small parameter e 
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one obtains that the equation is satisfied identically in 
the orders e° and e 1 . Next one has 

vd^(f>i — 2U\q\ — 2 cos 2 a q 2 = , . . 

vd ( (j) 2 - 2 sin 2 aq x - 2U 2 q 2 = ^ ' 

in the order e 2 and 

9 c 2 0i - vd cgi = 0, d 2 c <p 2 - vd c q 2 = (A2) 

I 



in the order e 3 . Integrating equations (|A2|) once and tak- 
ing into account l|16|). one obtains Eq. ((T7|l . Substituting 
link (H3 in (53J one verifies that the so obtained system 
is solvable subject to the condition 1|1U|I . what justifies 
the value of the soliton velocity as the group velocity of 
the sound waves. 

Equations of the order e 4 , where the link 117fl is used, 
read 



1 3/Ti 2 COS 2 CV COS 2 Cv 

vdtfu - 2U mi - 2cos 2 a g21 = d T ^ - -ftfa + ^(d c 0i) 2 + — (W(c^ 2 ) + ^^{d^ 2 f , 

V s V 2 $ 

vc\0 21 -2sin 2 agH - 2[/ 2 q 21 = T dp 2 - ^ 3 2 + ^ 2 -(9 c 2 ) 2 + ^_^(9 C ^)(5^ 2 ) + ^^(c^i) 2 ■ 



(A3) 



Finally one computes the equations of the order e (where 
the link (|17|l as well as the explicit value of the group 
velocity (|1(J[1 are taken into account and integration with 
respect to Q is performed) 

3 1 

f \ (A4) 
<9 c 2 i - vq 2 \ = -■^-(c'c'fe) 2 d T <j) 2 . 

The last system allows one to expresses qj\ through 
other dependent variables, and substitute the result in 



equations (|A3I) . In this way one obtains the inhomo- 
geneous linear algebraic system of equations for d^fiji, 
which determinant is computed to be zero. This leads 
to linear dependence of the equations which are solv- 
able only subject to the respective requirement imposed 
on their right hand sides. Differentiating the mentioned 
solvability condition with respect to £ and using one more 
time the link (|17|l in order to express the final result 
through the dependent variables qj only, one obtains cou- 
pled KdV equations (JTSJ. 
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